home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Technotools
/
Technotools (Chestnut CD-ROM)(1993).ISO
/
lang_c
/
cug188
/
translib.for
< prev
Wrap
Text File
|
1985-08-21
|
4KB
|
161 lines
/*
HEADER: ;
TITLE: C elementary transcendentals;
VERSION: 1.0;
DESCRIPTION║ááá Sourcσááá codσááá fo≥áá al∞áá standarΣááá ├ ì
transcendentals
Employ≤á ldexp(⌐á anΣá frexp(⌐á functions╗á iµ ì
suitablσá version≤á oµ thesσ arσá no⌠á provideΣ ì
b∙á ß giveε compiler¼á thσ version≤ provideΣ iε ì
sourcσá codσá wil∞ requirσá adaptatioεá t∩á thσ ì
doublσ floa⌠ format≤ oµ thσ compiler.
KEYWORDS: Transcendentals, library;
SYSTEM: CP/M-80, V3.1;
FILENAME: TRANS.C;
WARNINGS║áá frexp(⌐áá anΣá ldexp(⌐á arσáá implementatioε ì
dependent«á Thσá compile≥á employeΣá doe≤á no⌠ ì
suppor⌠ááá minu≤áá (-⌐áá unar∙áá operator≤áá iε ì
initialize≥á lists¼á whicΦ arσ requireΣ b∙á thσ ì
code.
AUTHORS: Tim Prince;
COMPILERS: MIX C, v2.0.1;
*/
program transt
real log
sum=0
do 1 i=1,1000
x=i*.01
1 sum=abs(tan(atan(x))/x-1)+sum
write(3,2)sum
2 format(' built-in tan(atan) error:',g15.7)
sum=0
do 3 i=1,1000
x=i*.01
3 sum=abs(exp(log(x))/x-1)+sum
write(3,4)sum
4 format(' built-in exp(log) error:',g15.7)
sum=0
do 5 i=1,1000
x=i*.01
5 sum=abs(xtan(xatan(x))/x-1)+sum
write(3,6)sum
6 format(' tan(atan) error:',g15.7)
sum=0
do 7 i=1,1000
x=i*.01
7 sum=abs(xexp(xlog(x))/x-1)+sum
write(3,8)sum
8 format(' exp(log) error:',g15.7)
end
function xtan(x)
xtan=xsin(x)/xcos(x)
return
end
function xcos(x)
xcos=xsin(1.57079633-x)
return
end
function xsin(x)
real table(5)
data pi/3.141592653589793238/,
& table/2.60198e-6,-.0001980746,.0083330258,
& -.16666657,.99999999/
c series to 28 bit accuracy
c nearest multiple of pi
pim=anint(x/pi)
c get remainder
xr=x-pi*pim
c sin(x-pim*pi) = (-1)**pim * sin(x)
hpim=scal2(pim,-1)
c if pim/2 not integer, pim is odd so change sign
if(aint(hpim).ne.hpim)xr=-xr
c evaluate Horner polynomial, odd terms
xsin=poly(xr*xr,5,table)*xr
return
end
function scal2(x,i)
c scal2=x*2**i
integer*1 ix(4),i
equivalence(xi,ix(1))
xi=x
ix(4)=i+ix(4)
scal2=xi
return
end
function poly(x,n,table)
real table(n)
poly=table(1)
do 1 i=2,n
1 poly=poly*x+table(i)
return
end
function xatan(x)
logical invert
real table(8)
data table/-.0046930,.0242506,-.0594827,.0991394,
& -.1401932,.1996969,-.3333199,.9999999/
xr=abs(x)
invert=xr.gt.1
if(invert)xr=1/xr
xr=poly(xr*xr,8,table)*xr
c tan(1/x)= tan(pi/2-x)
if(invert)xr=1.57079633-xr
xatan=sign(xr,x)
return
end
function xlog(x)
xlog=alog2(x)*.693147181
return
end
function alog2(x)
real table(3)
integer*1 ix(4),iexp
equivalence(xr,ix(1))
c polynomial to 24 bit precision
data table/.5957779,.9615884,2.8853901/
xr=x
c z'35' is leading 7 bits (minus hidden bit) of sqrt(.5)
c shift so xr close to 1.
c .true. value -1 on this system, bias z'80'
iexp=(ix(3).le.z'35')+ix(4)-z'80'
c subtract 1 first so no accuracy lost
xr=scal2(xr,-iexp)-1
xr=xr/(xr+2)
c polynomial in (xr-1)/(xr+1), .7 < xr < 1.4
alog2=poly(xr*xr,3,table)*xr+iexp
return
end
function xexp(x)
xexp=exp2(x*1.442695040888963407)
return
end
function exp2(x)
integer*1 expn
real table(7)
c 26 bits precision -.5 < xr < .5 for 2**x
c sinh (|x|<.7) use odd terms (x*1.442695)
data table/.000154,.00133904,.00961814,.05550358,
& .2402265,.69314718,1./
if(x.ge.-128)go to 1
exp2=0
return
1 xr=amin1(x,126.9999)
expn=anint(xr)
xr=xr-expn
exp2=scal2(poly(xr,7,table),expn)
return
end